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Abstract 

Biological  systems  usually  consist  of  large  number  of  components  and  involve  processes  at  a  variety  of 
spatial,  temporal  and  biological  scales.  Systems  biology  aims  to  understand  such  systems  by  integrating 
information  from  all  functional  levels  into  a  single  dynamic  model.  This  study  discusses  model  development 
and  the  use  of  Global  Sensitivity  Analysis  (GSA)  in  systems  biology  modeling  and  shows  how  the 
information  content  of  clinical  data  from  Short  Insulin  Tolerance  Test  (S1TT)  can  be  handled  by  optimal 
model-based  estimation  techniques.  The  goal  is  to  develop  dynamic  model  for  type  II  diabetes  and  estimate 
set  of  parameters  of  the  model  with  greater  accuracy  and  precision.  Based  on  the  SITT  data,  the  blood  glucose 
dynamic  model  was  developed  as  a  system  of  linear  differential  equations  with  constant  coefficients 
(parameters).  The  sensitivity  of  the  parameters  was  tested  using  a  novel  GSA  based  approach,  Derivative- 
Based  Global  Sensitivity  Measures  (DGSM).  The  proposed  approach  was  implemented  in  SensSB  (a  Matlab 
based  toolbox).  For  the  purpose  of  comparison,  the  sensitivity  of  the  model  was  also  tested  using  Sobol’s 
method  and  a  local  approach.  The  results  have  shown  that  the  model  is  less  sensitive  to  the  third  parameter 
Kn  and  fits  the  SITT  data  satisfactorily. 

Keywords:  Sensitivity  analysis,  Biological  systems,  SITT  data,  Blood  glucose  concentration  model,  Cramer- 
Rao  statistical  method 

1.  Introduction 

Sensitivity  analysis  (SA)  generally  is  the  study  of  how  the  uncertainty  in  the  output  of  a 
mathematical  model  or  system  can  be  apportioned  to  different  sources  of  uncertainty  in  its  inputs 
(Pannell  1997).  The  SA  methods  are  used  to  determine  the  influence  of  uncertainty  factors  in  the 
output  of  a  model.  They  also  indicates  which  of  the  uncertainty  factors  need  to  be  further  examined 
in  order  to  reduce  uncertainty  (Saltelli  and  Tarantola,  2006).  SA  algorithms  were  applied  in  physics 
(Pastorelli  et  al.,  2000),  chemistry  (Saltelli  and  Tarantola,  2006),  social  sciences,  and  many  others. 
However,  fewer  techniques  have  been  exploited  in  the  field  of  biomedical  engineering  (Hu  and  Shi, 
2010).  SA  methods  are  classified  into  local  methods  (LSA)  and  Global  methods  (GSA).  In  LSA, 
inputs  are  varied  one  at  a  time  by  a  small  amount  around  some  fixed  point  and  the  effect  of 
individual  perturbations  on  the  output  is  calculated.  GSA  on  the  other  hand,  allows  variation  of  all 
inputs  simultaneously  over  their  entire  input  space.  Specifically  GSA  uses  a  sampling-based 
approach,  and  the  effects  on  the  output  of  both  individual  inputs  and  interactions  between  inputs  are 
assessed  (Sumner  2010).  The  advantageous  characteristic  of  GSA  methods  is  that  they  are 
independent  of  model  characteristic  such  as  linearity  between  the  outputs  and  inputs  of  the  model 
(Buis  et  al.,  2010). 

SA  was  used  in  mathematical  modeling  in  many  areas  including  biological  model  fonnulations 
(Hetherington  et  al.,  2006).  The  most  common  example  of  SA  in  biological  system  is  the  use  of 
metabolic  control  analysis  (MCA)  in  the  study  of  metabolism  (Sumner  2010).  However,  it  was  less 
commonly  used  in  other  areas  of  biology  such  as  cellular  signaling  (Hu  and  Yuan,  2006). 
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Model  identification  system  describe  a  dynamic  behavior  of  a  system  or  process  in  either  the  time  or 
frequency  domain.  Most  of  the  models  used  in  blood  glucose  concentration  identification  are  Auto- 
Regressive  models  with  exogenous  inputs  (ARX)  to  represent  unknown  dynamic  of  blood  glucose 
concentration  (Isuru  et  al.,  2015).  Recently,  adaptive  system  identification  techniques  have  been 
used  for  prediction  of  blood  glucose  concentration  to  minimize  hypoglycemia  (Eren-Oruklua  et  al, 
2012).  In  this  paper,  a  dynamic  model  was  identified  and  sensitivity  analysis  technique  was  applied 
on  blood  glucose  concentration  model  for  type  II  diabetes  with  the  view  to  estimating  a  set  of 
parameters  of  the  model  with  reasonable  degree  of  accuracy  and  precision. 

2.  Material  and  Methodology 

The  Short  insulin  Tolerance  Test  (SITT)  was  conducted  at  the  General  Out-Patients  Department  (GOPD) 
of  University  of  Lagos  Teaching  Hospital  by  Fasanmade  (1987),  adhering  strictly  to  regulation  for 
use  of  human  subjects  as  reported  by  Gutti  et  al.  (2010),  Unlike  in  Gutti  et  al.  (2010),  here,  the  data 
was  classified  into  parameters  and  control  variables  and  stored  in  matrix  G  with  the  amount  of 
insulin  injected  as  input  variable  (manipulated  variable).  Also,  the  data  were  rearranged  and 
reshaped  into  a  single  column  such  that  it  can  represent  the  output  of  the  model. 

The  procedure  followed  in  model  development  and  sensitivity  analysis  on  blood  glucose 
concentration  model  is  shown  in  Figure  1. 


SITT  DATA 


Figure  1:  Procedure  for  Model  development  and  sensitivity  analysis. 

2.1  Problem  Generation 

The  interaction  of  blood  glucose  concentration  in  human  body  can  be  formulated  based  on  the  law  of 
conservation  of  mass,  as  follows: 

Rate  change  =  in  -  out  (1) 

From  Equation  (1)  above,  the  rate  change  of  glucose  concentration  can  be  written  as 

Rate  change  of  glucose  concentration  =  Glucose  in  -  Glucose  out  (2) 
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The  rate  change  of  glucose  concentration  is  determined  by  insulin-independent  glucose  uptake  in 
which  the  SITT  was  conducted  on  patients  with  non-insulin  dependent  diabetes  mellitus  within  the 
age  range  of  25-60  years  as  reported  in  Gutti  et  al.  (2010). 

Let  Gx  denote  the  amount  of  glucose  concentration  in  the  body  and  G2  denotes  the  blood  glucose 
concentration  at  time  t  >  0  .  From  the  information  available,  assuming  the  reaction  kinetics  follows 
first  order,  Equation  (2)  can  be  mathematically  written  as: 


dGx 

dt 


~k\Gx 


(3) 


d-^  =  KGx-k2(G2-GQ)  (4) 

dt 


where:  Gx  (mg/dL)  is  the  glucose  concentration  in  the  body,  G2  (mg/dL)  is  the  plasma  blood 
glucose  concentration,  G0  is  the  initial  blood  glucose  concentration  during  the  fasting  period 
(baseline  value)  at  t  —  0  .  The  initial  condition  are  G,  (o)  =  G2  (o)  =  G0 ,  the  rate  constant  K]  ( min  1 )  is 
the  rate  of  glucose  absorption  and  K2  (  min-1  )  is  the  rate  of  disappearance  of  glucose  and  t  is  the 
time  taken  during  the  experiment  (min).  Equations  (3)  and  (4)  represent  the  glucose-insulin  system 
following  the  short  insulin  tolerance  test  (SITT).  Note  that  in  Equation  (3),  the  short  insulin 
tolerance  test  was  done  during  fasting  and  there  was  no  flow  of  glucose  into  the  body. 

In  order  to  determine  the  blood  glucose  concentration  level  G2  (t)  explicitly,  we  must  first  solve  the 
differential  Equation  (3)  and  obtain  an  expression  for  Gj(t)  .  The  result  of  the  separation  of  the 
variables  in  Equation  (3)  is; 

-  Kxdt  =>  ln|G[ |  =  -Kxt  +  d 

Gx{t)=G0eK't  (5) 


Note  that  Equation  (3)  is  a  first  order  linear  differential  equation  of  the  form: 


dM+p(,)y(t)=q(,) 


It  has  been  verified  in  Ganesh  and  Balasubramanian,  (2009)  that: 


-f  p(t)dt 

y  =  e  J 


|  q(t)e^  P^dt  dt  +  C 


(6) 

(V) 


Equation  (7)  represent  the  solution  to  equation  (6)  where  C  is  the  constant.  Substituting  (5)  into  (4) 
gives  a  linear  first  order  differential  equation: 

^dll  +  KiGi(,)  =  Ki(Gae-^)+K2G0  (8) 

Equation  (6)  above  has  the  following  form  of  solution: 
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G2(t) 


+  K2G„)e^^dt  +  C 


G2(t)-Ga-^-r(e~'c'-e-,‘-')+GQ  (9) 

For  SITT  data,  the  insulin  is  administered  intravenously  and  therefore  the  response  of  plasma 
glucose  will  be  very  fast,  making  K2  »>  Kl  att  >  0 .  Based  on  this,  the  above  equation  reduces  to: 

G2{t)  =  G0e-Kt+G0  (10) 

Since  the  level  of  glucose  changes  with  respect  to  the  type  of  meal  taken,  it  is  represented  by  a 
system  of  linear  differential  equation.  One  major  challenge  in  system  identification  is  to  estimate  the 
unknown  parameters  and  perform  an  SA  to  verify  regions  of  the  system  and  the  most  relevant  state 
variables  (Rodriguez -Fernandez  and  Banga,  2010).  Estimating  the  unknown  parameters  of  a 
mathematical  model  requires  the  input-output  data  and  the  class  of  model  (Jacobs  2015).  The 
parameters  are  chosen  or  guessed  so  that  the  output  of  the  model  is  the  best  match  with  respect  to  the 
experimental  data  (Jacobs  2015). 

2.2  Model  Implementation 

SensSB  software  toolbox  was  used  for  the  model  implementation.  The  model  was  implemented  in 
1.5  GHz  core  i3,  4.00GB  RAM  system.  The  study  was  begun  with  estimation  of  parameters  then 
sensitivity  analysis  with  local  method  then  followed  by  global  approaches  based  on  DGSM  and 
Sobol’s  method  respectively. 

The  SITT  data  was  collected  and  rearranged  into  a  single  column  such  that  it  represented  the  output 
of  the  model.  The  algebraic  model  Equation  (9)  describing  the  SITT  was  developed  and  presented  in 
section  2.1  above.  In  order  to  implement  the  model  in  SensSB  software,  it  is  necessary  to  formulate 
the  problem  so  that  it  can  accommodate  the  experimental  data.  The  model  equation  was  solved  with 
the  MATLAB  solver  ode  15s  and  local  sensitivities  were  calculated  with  SENS  SYS.  In  the  ordinary 
differential  equations  (ODEs)  initialization  stage,  the  necessary  information  to  compute  the  function 
was  provided,  the  SITT  data  had  one  state  variable  (plasma  glucose  concentration  level)  and  the 
initial  range  for  nonnoglycemic  subject  was  taken  as  70-80  mg/dL  according  to  Diabetes 
Association  of  Nigeria  (Chinenye,  2013).  The  data  had  10  sampling  points  at  [-1002468  10  12  14 
16]  measured  in  minutes.  The  system  of  ODEs  involves  3  parameters,  which  are  G0  ,  K]  ,  and  K2  . 

They  were  represented  by  pi,  p2,  and  p3  respectively.  The  initial  value  chosen  for  the  three 
parameters  are  [75  0.028  0.026],  lower  bound  [60  0.01  0.01]  and  upper  bound  [80  0.03  0.03].  The  3 
parameters  are  being  estimated  so  the  vector  [1:3]  equivalent  to  [1  2  3]  were  introduced  as  ‘Index  of 
the  parameters  to  be  modified’,  please  mind  that  the  upper  and  lower  bounds  are  going  to  influence 
the  value  of  the  resulting  global  sensitivity  analysis. 

3.  Results  and  Discussion 

The  results  of  the  model  identification,  SA  and  Parameters  estimation  are  presented  and  discussed  in 
this  section. 
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3.1  Dynamic  model 

The  results  obtained  for  the  dynamic  model  described  in  Equations  (3)  and  (4)  was  fitted  on  SITT 
data  and  shows  an  appreciable  consistence  with  the  level  of  blood  glucose  concentration  measured 
during  SITT  experiment.  The  dynamic  model  results  obtained  was  compared  with  Estela,  (2011) 
results  which  indicates  a  better  fit  on  the  SITT  data.  Figures  2  and  3  shows  a  comparison  between 
the  model  predicted  values  and  the  SITT  data  reported  by  Gutti  et  al.  (2010),  for  a  control  and 
diabetic  subject  corresponding  to  the  blood  glucose  concentration. 


Figure  2:  Fitted  SITT  data  and  the  dynamic  model  for  a  control  subject  in  SensSB. 


Figure  3:  Fitted  SITT  data  and  the  dynamic  model  for  diabetic  subject  in  SensSB 

It  can  be  observed  that  the  estimated  parameters  allow  reproducing  almost  exactly  the  experimental 
data  while  all  the  local  solvers  that  were  tried  with  this  initial  point  failed  to  converge,  or  converged 
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to  bad  local  solutions.  Different  global  solvers  such  as  stochastic  ranking  evolutionary  search  (SRES) 
or  differential  equation  (DE)  also  failed  or  converged  in  a  much  larger  computational  time. 

3.2  Sensitivity  analysis  of  the  model 

Here,  the  results  for  the  sensitivity  analysis  of  the  model  based  on  local,  Sobol  and  DGSM  method 
are  presented.  These  are  shown  in  Tables  1,  2  and  3,  while  their  respective  profile  are  shown  in 
Figures  4,  5  and  6  respectively.  The  parameter  ranking  based  on  local,  Sobol’s  and  DGSM  method 
are  shown  Figures  7,  8  and  9  respectively. 


Table  1:  Local  sensitivity  analysis  results 

Ranking 

Parameter 

msqr 

Mabs 

Mean 

Min 

Max 

abs_sens 

1 

P2 

1.16e-001 

4.02e-002 

-4.02e-002 

5.58e-012 

-4.65e-001 

8.10e+001 

2 

PI 

8.78e-002 

4.29e-002 

-4.85e-003 

2.08e-001 

-2.88e-001 

1.59e+002 

3 

P3 

0.00e+000 

0.00e+000 

0.00e+000 

0.00e+000 

0.00e+000 

0.00e+000 

Table  2:  Global  sensitivity  analysis  results  based  on  Sobol  method 


Parameter 

sij 

SIJAT 

PI 

2.496e-001 

8.293e-001 

P2 

7.504e-001 

3.378e-001 

P3 

0.000e+000 

0.000e+000 

Table  3:  Global  sensitivity  analysis  results  based  on  DGSM  method 

Parameter 

abs_MJA* 

abs_sigmajA* 

rel  MjA* 

rel  sigmajA* 

PI 

6.198e-001 

6.160e-001 

6.501e-001 

6.219e-001 

P2 

3.802e-001 

3.840e-001 

3.499e-001 

3.781e-001 

P3 

0.000e+000 

0.000e+000 

0.000e+000 

0.000e+000 

where:  msqr  is  the  parameter  ranking  based  on  the  squared  root  of  the  square  relative  sensitivities, 
mabs  is  the  parameter  ranking  based  on  the  mean  of  the  absolute  values  of  the  relative  sensitivities, 
mean  represent  the  parameter  ranking  based  on  the  mean  of  the  relative  sensitivities,  min  is  the 
ranking  based  on  the  minimum  value  of  the  relative  sensitivities,  max  is  the  ranking  based  on  the 
maximum  values  of  the  relative  sensitivities  and  abs  sens  is  the  parameter  ranking  based  on  the 
mean  of  the  absolute  sensitivities. 
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Figure  5:  Sobol’  profile  for  dynamic  model. 
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parameter 


Figure  7:  Parameter  ranking  based  on  local  method 


Figure  8:  Parameter  ranking  based  on  Sobol’  Figure  9:  Parameter  ranking  based  on  DGSM 
method  method 

As  can  be  seen  from  the  Figures  7,  8  and  9,  parameter  3  shows  a  low  relative  sensitivity,  which 
indicate  the  model  is  less  sensitive  to  the  third  parameter.  This  is  observable  in  both  sensitivity 
analysis  methods.  Subsequently,  local  and  global  sensitivity  analysis  were  conducted  to  determine 
the  most  sensitive  parameter  of  the  model  identified  in  Section  2.1.  It  shows  that  the  method  can  be 
easily  implemented  in  MATLAB.  The  results  of  the  sensitivity  analysis  further  show  that  DGSM 
computations  are  reliable  to  a  reasonable  degree  of  accuracy,  while  Sobol’s  method  ranked  the 
parameters  in  order  of  importance. 

3.3  Parameter  estimation  of  the  model 

The  results  of  the  parameter  estimation  on  local,  Sobol’  and  DGSM  method  are  presented  in  Tables 
4,  5  and  6  respectively. 

Table  4:  Parameter  estimation  results  based  on  local  method _ 

Optimal  objective  function  value  (least  square)  =  1.460e+000 
CPU  time  required  for  the  estimation  =  3.24e+002sec 

Parameter  Optimal  value  Confidence  interval  (Cramer-Rao) 

PI  1.500e-002  +-2. 689-004 
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P2 

2.260e-002 

+-3.537e-004 

P3 

7.590e+001 

00 

Table  5:  Parameter  estimation  results  based  on  Sobol’  method 


Optimal  objective  function  value  (least  square)  =  1.216e+000 
CPU  time  required  for  the  estimation  =  3.00e+002sec 


Parameter 

Optimal  value 

Confidence  interval  (Cramer-Rao) 

PI 

1.607e-002 

+-5.190e-004 

P2 

3.000e-002 

+-7.232e-004 

P3 

6.970e+001 

00 

Table  6:  Parameter  estimation  results  based  on  DGSM  method 

Optimal  objective  function  value  (least  square)  =  1.423e+000 

CPU  time  required  for  the  estimation  =  3.01e+002sec 

Parameter 

Optimal  value 

Confidence  interval  (Cramer-Rao) 

PI 

1.620e-002 

+.4.931-004 

P2 

3.000e-002 

+-6.859e-004 

P3 

6.171  e+00 1 

OO 

The  parameter  estimated  based  on  DGSM  differ  from  those  obtained  using  local  and  Sobol’s  method 
because  DGSM  depend  on  nominal  values  of  the  parameter.  Estimated  parameter  values  varies  with 
different  subjects  and  contradicts  earlier  reports  by  Sumner,  (2010)  who  uses  the  principal 
components  analysis  to  measure  the  importance  of  a  parameter  on  the  entire  model  output.  The 
sensitivity  analysis  results  presented  in  this  work  are  based  on  the  dynamic  model  of  blood  glucose 
concentration  which  is  the  goal  of  this  study.  The  LSA  and  GSA  based  method  are  assessed  on  the 
basis  of  the  length  of  the  confidence  intervals  computed  for  the  parameters  with  Cramer-Rao 
statistical  method  (Rodriguez -Fernandez  and  Banga,  2010)  as  shown  in  Tables  4,  5  and  6  in  which 
the  DGSM  outperfonn  both  Sobol’s  and  local  methods. 

4  Conclusion 

The  prevalence  of  type  II  diabetes  is  increasing  in  Nigeria  and  globally  and  most  of  the  mathematical 
models  developed  were  devoted  to  the  dynamic  of  glucose-insulin  system  based  on  Intravenous 
glucose  tolerance  test  (IVGTT),  oral  glucose  tolerance  test  (OGTT)  and  frequently  sampled 
intravenous  glucose  tolerance  test  (FSIVGTT)  data.  This  study  proposed  Short  Insulin  Tolerance 
Test  (SITT)  based  model  via  optimal  design  of  experiment.  Our  approach  incorporates  important 
steps  required  in  model  development  process  such  as  global  sensitivity  analysis,  Pseudo-global 
identifiability  analysis,  Optimal  experimental  design  based  on  global  sensitivities,  robust  parameter 
estimation,  etc.  Techniques  of  global  sensitivity  analysis  appear  more  appealing  due  to  the  fact  that 
they  investigate  the  entire  range  of  parametric  uncertainty.  The  parameters  for  the  model  developed 
were  estimated  in  the  least  square  sense  with  SITT  data  fitting  the  model  satisfactorily.  The  result  of 
local  sensitivity  was  compared  with  Sobol’  and  DGSM  and  the  values  are  assessed  based  on  the 
length  of  confidence  intervals  with  Cramer-Rao  statistical  replication. 
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